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FOREWORD 


This  reporx  extends  a previously- 
developed  (Ref.  1)  grav'imetric  estimation 
algorithm  by  removing  a data  length  restric- 
tion associated  with  the  earlier  algoritim. 

The  original  method,  which  could  efficiently 
process  very  large  amounts  of  data,  required 
that  measurements  be  available  over  an 
interval  much  longer  than  the  correlation 
distances*  associated  with  the  gravimetric 
quantities  being  measured.  This  restriction 
could  be  significant  in  problems  where  long 
correlation  distances  occur  such  as  processing 
geoid  height  data. 


In  the  current  study  modern  signal 
processing  theory  involving  discrete  Fourier 
transforms 
estimation 
is  an  effi 

cable  to  processing  large  quantities  of 
The  new  algorithm  is  called  GEOF.riST , an  acr_pjiym 
for  Geodetic  Fast  Estimation.  {>■[ 


is  used  to  solve  the  finite  length 
problem.  The  result  of  tne  effort 
:ient  algorithm  particularly  appli- 

data. 
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1.  INTRODUCTIOy 

1 . 1 BACKGROUND 

Continuing  developments  in  gravimetric  sensor  tech- 
nology are  resulting  in  the  availability  of  large  amounts  of 
gravity  field  data.  In  addition  to  greater  densif ication  of 
surface  gravity  measure.ments , large  holdings  of  GEOS-3 
satellite  alt imetric  data  are  now  available.  Measurements  of 
the  ocean  gecid  will  increase  in  abundance  with  the  orbiting 
of  SEASAT-A.  The  advent  of  moving-base  gradiometry  and 
utilization  of  satellite-to-satellite  tracking  measurements 
will  add  even  more  to  the  applicable  data  base. 

The  gravity  field  quantities  which  must  be  extracted 
from  this  heterogeneous  data  base  are  equally  varied.  Values 
of  the  gravity  disturbance  vector  are  needed  to  compensate 
terrestrial  inertial  navigation  and  guidance  systems.  Gecid 
undulations  are  needed  for  charting  ocean  currents  and  for 
determining  elevation  in  geodetic  surveys.  Estimates  of 
spatial  derivatives  of  gravity  are  required  for  gravity 
gradiometer  testing  and  e\*aluation.  A schematic  which  illus- 
trates exa.T.ples  of  gravity  field  measurements  of  interest  and 
lists  some  of  the  gravimetric  quantities  to  be  estimated  from 
such  measurements  is  presented  in  Fig.  1.1-1.  These  and  other 
geophysical  applications  have  the  requirement  for  rapid,  effi- 
cient grari.metric  data  processing. 

Because  of  interrelations  among  all  of  the  gravity 
field  quantities,  and  because  the  value  of  the  field  at  each 
point  is  signif icantly  correlated  with  the  gravity  field  at 


More  recently,  statistical  techniques  have  been 
applied  to  optimally  estimate  gravity  field  quantities  (Hef.  3) 
and,  in  many  instances,  have  provided  considerable  advantage 
over  classical  methods  (Refs.  4 and  7).  One  feature  which 
makes  the  statistical  approaches  attractive  is  "built-in 
error  analysis"  which  allows  the  impact  of  data  deficiencies 
such  as  noise,  discretization  and  finite  extent  to  be  easily 
studied.  For  some  proble.ms  involving  gravity  effects  along 
a single  track  (such  as  .might  be  traversed  by  an  inert ially 
navigated  terrestrial  vehicle),  recursive  methods  such  as 
Kalman  filtering  have  proved  useful  for  processing  large 
amounts  of  data  (Ref.  5).  However,  the  recursive  approaches 
rely  on  state  space  formulations  which  are  not  easily  applied 
to  many  gravimetric  problems,  particularly  those  which  involve 
two  or  three  dimensions.  As  a result,  for  most  gravity 
quantity  estimation  problems,  "classical"  least  squares 
techniques  must  be  applied.  Because  such  methods  involve 
matrix  inverses,  excessive  computational  burdens  result  when 
a large  number  of  measurement  points  is  used  (say  1000  or 
more).  Of  course  a least  squares  estimation  problem  can 
always  be  truncated  to  reduce  it  to  manageable  size  (e.g., 

Ref.  6),  but  the  error  involved  in  such  approximations  can 
be  unacceptable.  Ivaluation  of  the  error  may  also  be  dif- 
ficult . 

Another  approach  to  least  squares  gravimetric  data 
processing  is  developed  in  Ref.  1.  The  concept  exploits  the 
special  structure  of  the  gravity  model  covariance  matrices 
instead  of  approximating  either  the  data  or  the  integral  rela- 
tions a.mong  gravimetric  quantities.  The  formulation  takes 
advantage  of  the  computational  efficiency  of  fast  Fourier 
transforms  and  results  in  an  estimation  algorithm  which  is 
suitable  for  processing  gravimetric  data  sets  containing  very 
large  nu.mbers  of  measurements  (100,000  measurements  provide 
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no  dif  f iculr  y ) . However,  the  algoritar.  developed  in  Ref.  1 
is  subject  to  moderate  levels  of  error  when  the  data  interval 
is  not  significantly  larger  than  the  correlation  distance  of 
the  gravity  quantities  being  measured. 


1.2  PURPOSE  .4XD  SCOPE  OF  THIS  REPORT 

In  this  report,  further  development  of  the  frequency 
domain  gravimetric  estimation  method  developed  in  Ref.  1 is 
presented.  In  particular,  the  algorithm  is  extended  to  elimi- 


nate  errors  cue  to  an  ancrc.xtmati 


optimal  and  still  retains  the  necessary  computational  soeed 
recuired  to  solve  very  high  order  gravimetric  estimation 
problem.s . Of  course,  as  with  any  estimation  technique,  errors 
arising  from  finite,  noisy  data  of  lim.ited  extent  still  cause 


esti.m.ation  er: 


•s.  The  merit  cf  the  estim.ation  technique 


descrioec  in  this  report  is  that  no  computational  errors  are 
added  in  the  course  of  processing  the  data. 


GEOFAST  fully  and  explicitly  takes  into  account  the 
finite  range  of  data  available  from,  most  gravimetric  surveys 
without  limiting  the  generality  cf  problems  to  which  the 
algorithm  is  applicable.  Thus,  the  augmented  algorithm  is 
suited  to  m.ultisensor  gravimetric  estimation  problem.s  (as  well 
as  those  involving  a single  measured  quantity)  in  exactly  the 
same  fashion  as  the  earlier  approach.  Specifically,  suitable 


'Results  in  niriim.um  variance  of  the  estimation  error  when 
correct  statistical  gravimetric  models  are  utilized. 

-The  term  gravimetric  is  used  in  the  broad  sense  of  relating 
to  techniques  for  measuring  or  estimating  the  gravitational 
disturbance  potential  or  any  of  its  derivatives. 


r 


data  inputs  include  anj’  gravimetric  quantity  available  at 
gridpoints  on  or  above  the  earth's  surface.  The  purpose  of 
this  report  is  to  describe  enhancements  to  the  frequency 
domain  estimation  algorithm  of  Ref.  1 which  result  in  an 
optimal  procedure  for  estimating  gravimetric  quantities  from 
multisensor  data. 

In  Chapter  2,  the  mathematical  development  of  the 
algorithm  is  presented  in  some  detail.  Chapter  3 sum.marizes 
the  enhanced  algorithm,  and  indicates  the  logical  ne.xt  steps 
for  its  exploitation.  Recommendations  for  useful  tests  of 
GZOFAST  are  also  presented. 

1 

i 

l. 3  TECHXIC.i.L  .4PPR0.4CH 

.Minim.u.T  variance  estimation  as  applied  to  gravimetric 
quantities  is  described  in  Ref.  1.  TTithout  repeating  that 
discussion,  it  is  appropriate  to  review  the  rationale  for 
utilizing  the  frequency  domain  formulation,  namely  that  the 
burden  of  manipulating  large  covariance  m.atrioes  associated 
with  the  esti.mation  process  is  eased.  The  discrete  Fourier 
transformation  takes  advantage  of  the  special  structure  of 
covariance  matrices  involved  in  the  estim.ation  problemi.  The 
matrices  that  require  inversion  in  the  discrete  Fourier 
transform  formulation  are  diagonal  or  band-diagonal.  In  the 
case  of  finite  data  lengths,  the  number  of  superdiagonal 
bands  (and  hence  the  computation  effort  required  to  effect  a 

m. atrix  inverse),  is  proportional  to  the  sidelobe  energy  of 

’ the  finite  length  data  spectrum..  To  minim.ize  this  sidelobe 

energy,  and  thereby  m.inim.ize  the  num.ber  of  superdiagonal  bands 
in  the  spectral  density  matrix,  a data  "window"  is  used.  The 
GEOFAST  algorithm,  consists  of  transforming  the  "windowed"  data 
into  the  discrete  frequency  dom.ain  and  there  solving  the 
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esximation  equations  using  correspondingly  windowed  and  trans 
fortred  covariance  functions.  The  actual  solution  of  the 
banded-diagonal  set  of  equations  in  the  frequency  domain  is 
accomplished  by  the  Cholesky  decomposition  technique.  An 
overview  of  the  complete  approach  is  depicted  in  Fig.  1.3-1. 


niotn 
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T 


Figure  1.3-1  GEOFAST  Algorithm 
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2.  FAST  COMPUTATIONAL  METHODS  FOR  ZSTI.\LATION 

2.1  COMPUTATIONAL  ELEMENTS  OF  GEOFAST 

The  various  statistical  procedures  for  gravimetric 
estimation  have  several  mathematical  operations  in  common. 

The  most  important  example  of  such  an  operation  is  the  solu- 
tion of  a system  of  equations  involving  the  data  covariance 
matrix.  k second  i.mpcrtar-t  operation  is  the  vector-matrix 
multiplication  typicallv  carried  out  after  solving  the 
covariance  equations.  These  straightforward  algebraic  opera- 
tions are  significant  because  of  the  large  amounts  of  computer 
time  they  can  require.  For  many  large  gravity  quantity 
estimation  problems,  a conventional  implementation  of  the 
estimation  equations  leads  to  prohibitive  computer  processing 
requirements.  The  computational  load  occurs  primarily  because 
standard  algebraic  aigorith.ms  do  net  take  account  of  the 
special  mathematical  structure  found  in  gravity  quantity 
esti.mation  problems.  This  chapter  is  concerned  with  fast 
computational  .methods  for  carrying  out  the  fundamental  mathe- 
matical operations  in  gravity  quantity  estimation.  In  par- 
ticular a technique  which  exploits  the  special  structure  of 
the  gravity  estimation  covariance  matrices  is  developed. 

The  class  of  statistical  procedures  for  gravity  quan- 
tity estimation  to  which  the  techniques  of  this  report  are  most 
applicable,  has  been  called  the  "method  of  least  squares 
collocation"  (Ref.  3).  This  method  utilizes  minimum  variance 
estimation  for  a model  of  the  form 

z = Kx  V .2.1-1 


5 


• BLUE  (best  least-squares  unbiased  estir.ation) 


I If  X is  assumed  to  be  a stochastic  vector  and  H = I, 

^ where  I is  the  unit  matrix,  then  an  estimation  procedure  is 

I obtained  which  is  referred  to  here  as  im.plicit  minimum,  variance 


estimation 

(IMVE).  The 

terminology  em.phasizes  that 

the  rela- 

tion  between  x and  z is 

given  im.plicitly  through  a 

crcss- 

covariance 

miatrix  C 

If  X is  assumed  to  be  a det 

erm.inist  ic 

xz 

vector  and 

H defines  an 

explicit  (linearized)  funct 

ional  re- 

lation  between  x and  z,  the  procedure  is  termed  best 
_least-SGuares  unbiased  estim.ation  (BLUE).  This  is  sim.ply 
another  name  for  standard  least-squares  parameter  estimation. 
Reference  3 shows  that  if  x is  partly  stochastic  and  partly 
deterministic,  the  estimiation  procedure  decouples  into  the 
two  cases  just  defined.  A third  form  of  the  method  results  if 
X is  assumed  to  be  stochastic,  but  related  explicitly  to  z 
through  the  mieasurement  m^atrix  K.  This  case,  which  is  desig- 
nated explicit  minimium,  variance  _estim>ation  (EMVE^,  has  been 
shown  (Ref.  6),  to  be  equivalent  to  the  first  case  if  x is 
given  an  a nriori  covariance  miatrix  C and  the  distribution 

XX 

of  z is  defined  through  Eq . 2.1-1.  This  case  corresponds  to  a 
standard  Flalm.an  filter  implementation  (Ref.  19). 


A 
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2.1.1  IM\T 


The  mathematical  operations  of  IMl^are  representative 
of  those  found  in  gravity  quantity  estimation  in  general.  In 
this  formulation  (Ref.  1),  the  estimate  vector,  x,  is  found 
from  the  data  vector,  z,  via 


X 


(2.1-2) 


where  C and  C are  covariance  matrices  obtained  from  a 
zz  xz 


self-consisf-ent  statistical  gravity  model.  If  the  di.mensicns 
of  X and  z are  m and  n,  respectively,  then  and  C^._  are 

respectively  n^n  and  m>:n  matrices.  The  dimensionality  of 
these  'two  matrices,  and  their  analogs  in  other  formulations, 
is  the  source  of  the  computational  burden  in  estimation  prob- 
lems: for  large  n and  m,  the  implementation  of  Eg.  2.1-2 

requires  very  large  amounts  of  computer  time  and  is  subject  to 


severe  nrocessing  errors.  The  calculation  of  the  vector  C 


- j. 


zz  — 

.s  the  dominant  cart  of  this  calculation.  Efficient  eeneral 


purpose  algorithms  for  this  operation  require,  at  the  minimum. 


a number  of  numerical  operations  (multiplications  anc  ad 


r"  -ions; 


3 , 


which  is  approxi-mately  n /3 . Table  2.1-1  illustrates  the  rapid 


growth  of  processing  time  wizh  increasing  data  size  that  t 


cubic  law  produces.  The  operation  which  completes  Eq . 2.1-2 
is,  of  course,  the  multiplication  of  the  vector  z by  the 

cross-covariance  matrix  C . Straightforward  implementation 
of  this  operation  requires  mn  multiplications  and  m(n-l)  addi- 


tions. vrhen  the  number  of  quantities  to  be  estim.ated,  m,  is 
small  this  operation  is  relatively  inexpensive  and  may  be 
carried  out  by  standard  algorithms.  'When  m approaches  or  sur- 
passes the  order  of  n,  a special  purpose  algorithm  for  co- 
variance  matrix  multiplication  can  significantly  reduce  the 
cost  of  this  operation.  The  first  diagram  in  Fig.  2.1-1 
illustrates  the  decomposition  of  Eq . 2.1-2  into  the  operations 
of  covariance  equation  solution  and  covariance  matrix  multi- 


A 
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TABLE  2.1-1 


CO.MPUTATIONAL  COST  OF  .MINIMUM  VARIANCE  ESTIMATION 


OPTIMISTIC 

DATA 

estimate  of 

POINTS 

PROCESSING  TIME 

SINGLE  TRACK 

1000 

15  min 

5 - TRACK 

5000 

30  hr 

10  - TRACK 

10000 

10  days 

2.1.2  -EMIT  and  BLUE 

The  operations  associated  with  the  e.xplicit  .teasu 
ment  equation  formulation  of  EinT  and  with  parametric  BLUE 
are  similar  to  those  for  I.Ml’T . From  Refs.  3 and  6,  the 
expressions  analogous  to  Eq.  2.1-2  are: 


(EM^-E) 


(BLUE)  X = (H'C^.*  H)~^  h'c";  z 


(2, 


where  H is  an  n=<m  measurement  matrix,  C^,^.  is  the  n>^n  measu 
ment  noise  covariance  matrix,  and  C is  the  m^m  a nriori 
covariance  matrix  for  x.  The  superscript,  T,  indicates  ma 
transpose. 


Equations  2.1-3  and  2.1-4  may  both  be  convenient! 


decomposed 

into  three  sequen 

tial 

operat i 

ons 

• 

Computation  of 

a = 

- 

• 

.Multiplication 

to 

form  b = 

u' 

• 

Computation  of 
X = (H'c‘^  K ^ 

/V 

X = 

C"^ 

XX 

'-1  K 

H)' 
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MEASUngMFNIS 


IMPLICIT  minimum  variance  ESTIMATION  I IMVE  I 


EXPLICIT  minimum  variance  estimation  ( EMVE  1 


BEST  LEAST  SQUARES  UNBIASEC  ESTIMATION  ' BLUE  ' 


Figure  2.1-1  Block  Diagracrs  of  Gravi.merric 

EstitTiation  Algorithms 


Block  diagrams  of  the  EMVE  and  BLUE  processes  are  presented 
in  Fig.  2.1-1.  The  first  two  of  these  operations  are  subject 
to  the  same  considerations  and  the  sa.me  special  processing 
methods  as  are  the  operations  of  IMVE.  However,  for  many 
problem, s C . is  a diagonal  miatrix  (although  this  is  by  no 
means  always  the  case),  so  that  the  solution  of  the  corre- 
sponding covariance  equations  is  co.mputationaily  inexpensive. 
Vvhen , in  conjunction  with  diagonal  the  dimension  of  x is 

not  large,  both  the  EMVE  and  BLUE  method  m,ay  be  efficiently 
implemented  without  special  attention  to  cor,putat ionally  fast 
methods . 

Special  computational  methods  for  the  third  opera- 


tion associated  with  the  EMIT  or  BLUE  procedures  are  outside 
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the  scope  of  this  report,  vrhen  x is  not  of  large  dinensicr. , 

as  is  ordinarily  the  case  for  these  methods,  the  systems  of 

T _i  T -1  -1 

equations  involving  K C * H or  H C,,  . H + C are  not  large 
enough  to  require  special  numerical  attention.  However,  the 
efficient  com.putation  of  these  matrices,  which  involves  , 
is  made  possible  by  use  of  the  techniques  developed  in  this 
chapter  for  the  solution  of  covariance  equations.  Methods 
for  efficientlv  "inverting’’  K or  h”c~^  H when 

the  dimension  of  .x  is  large,  remain  a subject  for  future  work, 


2.1.3  Solution  Technicues 


The  remainder  cf  this  chapter  is  devoted  to  the 
description  and  analysis  cf  fast  computational  methods  for 
carrying  cut  the  two  principal  operations  which  have  been 
identified  above: 


Solution  of  large  systems  of  equations 
involving  or  C.„, 


Fast  multiolicat ion  cf  matrices  C 

X.: 

H when  the  dimension  of  x is  large 


or 


The  methods  discussed  in  this  chapter  are  those  appropriate 
for  gravity  quantity  estimation  problems.  Th'pically  such 
problems  entail  stationary  statistical  models  and  an  opera- 
tor relating  the  various  gravity  quantities  which  does  not 
change  under  spatial  translations. 

The  method  presented  in  the  pages  which  follow  for 


solution  of  stationary  covariance  equations  is 


1 . 3.n  c — s 


an  extension  of  the  classical  frequency  domain  Wiener  method 
The  basic  approach  of  generalized  Wiener  filtering  is 
described  in  Ref.  14  and  depends  for  its  success  on  the 
fast  Fourier  transform  (FFT)  technique  (Ref.  15).  The  use 
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of  classical  Wiener  filtering  techniques  in  conjunction  with 
gravity  data  processing  has  been  discussed  in  Ref.  1.  This 
report  develops  an  extended  computational  technique  which 
retains  the  speed  of  the  Wiener  method,  and,  at  the  sa-me 
time,  offers  the  high  degree  of  accuracy  required  for  gravity 
quantity  estimation. 


2.2  SPECTRAL  REPRESEXTATIO.V  OF  TOEPLITZ  MATRICES 

The  new  computation  methods  about  to  be  presented 
herein  are  based  on  so-called  Toeplitn  .matrices  (Ref.  Ic). 
Toepli.ta  matrices  consist  of  a class  of  n^n,  real-valued 
matrices  which  includes  the  covariance  matrices  from  sta- 
tionary random  processes.  More  formally,  a matrix  T is  of 
the  Toeplitz  type  if  there  is  a real-valued  function,  c., 
such  that  the  elements  of  T obey 


[T]  = C.  . , 0 < j , k < n-1 

j it-j 


(2.2-1) 


T is  an  autocovariance  or  a cross  covariance  matrix,  then 
is  the  corresponding  covariance  function.  The  m,atrix  T 
svmmetric  if  the  function  c.  is  even. 

The  relevance  of  Toeplitz  matrices  to  statistical 
estimation  has  already  been  hinted  at  in  the  definition  above. 
The  cova.-laace  matrices  (e.g.,  C,,^.  wblch 

appear  in  the  various  expressions  for  IMt’::.,  E’Ev E , and  BLUE 
are  of  the  Toeplitz  type  when  the  underlying  statistical 
models  are  stationary.  In  addition,  for  problems  involving 
the  esti.mation  of  gravity  quantities,  the  measurement  matrix 
H is  either  a Toeplitz  matrix  or  a submatri.x  fro.m  a Toeplitz 


'£ 

is 
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matrix.  This  follows  from  the  fact  that  the  geopotential 


operators  connecting  the  various  gravity  quantities  such  as 
the  Stokes,  Vening  Meinsez,  and  upward  continuation  opera- 
tors, are  invariant  under  spatial  shifts  (see,  for  instance 
Ref.  21).  Toeplitz  matrices  result  when  operators  with  this 
property  are  reduced  to  finite  dimensions,  typically  by 
discretization  and  truncation.  As  a consequence  of  sta- 
tionary gravity  models  and  shift-invariant  measurement 
operators,  essentially  all  large  matrices  appearing  in  the 
gravity  quantity  estimation  problem  are  of  the  Toeplitz  type. 
The  study  of  fast  computational  methods  for  gravity  quantity 
estimation  is  largely  the  study  of  efficient  .methods  tc 
handle  equations  defined  by  Toeplitz  matrices. 

9 

2.2.1  Fourier  Transformation  of  Circulant  Matrices 

Ciroulant  matrices  are  a subset  of  Toeplitz  m.atrices 
with  an  especially  si.mpie  form  under  Fourier  transf orm.ation . 

A Toeplitz  m.atrix  is  a ciroulant  (Ref.  S),  if  the  function 
defined  by  Ec . 2.2-1  obeys 


Matrices  with  this  property  m.ay  be  identified  by  the  fact 

that  each  row  is  equal  to  the  row  preceding  it  shifted  one 

element  to  the  right,  with  the  last  element  ’’wrapped  around" 

to  the  first  place.  For  example,  a circulant  matrix  with 

the  first  row  (c>,,  c, o ) necessarilv  has  its  second 

U 1 n-l 

row  given  bv  (c  , . c-,.  6-  , . . . , z r,)-  While  circulant 

^ • n-l'  O'  1 n-2 

m.atrices  themselves  do  not  appear  frequently,  their  properties 
are  useful  in  dealing  with  the  more  general  Toeplitz  m.atrices. 
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Circulant  matrices  have  the  special  property  that 
they  are  diagonalized  by  the  discrete  Fouriei  transformation 
(DFT).  The  n dimensional  DFT  matrix  F is  the  symmetric 

comp  lex- valued  matrix  whose  elements  are 


[P] 


jk 


~ exp 
/n 


2vi.-ik 

n 


0 <_  j , k t 


(2.2-3) 


and  its  inverse  is  easily  seen  to  be  F . The  DFT  of  an 
n-vector  of  complex  numbers,  x,  is  defined  by  the  relations 
(Ref.  9) 


X = F x,  X = P~X  (2.2-4) 

The  matrix  similarity  transformation  corresponding  to  Ec.  2.2-4 
is  given  by 


C = FCF 


(2.2-5) 


where  C is  an  n^m  matrix  of  complex  numbers  and 
representation  in  the  Fourier  transform  domain, 
ness  of  the  circulant  matrix  definition  depends 
following  result. 


C is  its 

Ui 

The  useful- 
largely  on  the 


The  Fourier  transformation,  Eq . 2.2-5,  of  any  circu- 
lant matrix  C is  the  diagonal  matrix 


= >'n  diag(OQ,  . . . , o^_^) 


(2.2-6) 


*The  symbol  i denotes  y -1 , and  the  superscript  (t)  denotes  the 
complex  conjugate  transpose. 
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T 


1 


where  £ is  the  discrete  power  spectrum 


4 = 


(2.2-7) 


and  4 = (o«,  . . . , 0 , corresponds  to  the  first  row  of  C. 

w li  “ X 

Moreover,  any  diagonal  matrix  is  the  Fourier  transform  of  a 
circulant  matrix,  defined  by  Eqs.  2.2-6  and  2.2-7.  If  the 
matrix  C is  symmetric  the  inverse  transform  F may  be  replaced 
by  F in  Eq.  2.2-7.  A proof  of  this  result  may  be  found  in 
Hef.  9.  This  theorem  is  also  equivalent  to  the  well  knowr 
relationship'  between  discrete  convolutions  and  discrete 
Fourier  transforms  proved  in  Ref.  10. 

2.2.2  Fast  Toeplitz  Matrix  Multiplication 

A computationally  fast  algorithm  for  multiplying  a 
vector  by  a Toeplitz  matrix  may  be  obtained  as  a consequence 
of  the  circulant  matrix  result  presented  in  Section  2.2.1. 

This  method,  together  with  the  fast  procedure  for  solving 
linear  equations  described  in  Section  2.3,  are  the  major 
components  of  an  efficient  gravity*  quantity  estimation  algo- 
rithm. Together,  the  two  algorithms  provide  a gravity  data 
processing  method  with  computer  time  requirements  proportional 
to  n log  n,  where  n is  the  number  of  data  points. 

In  gravity  quantity  estimation  the  need  for  an  opera- 
tion equivalent  to  multiplying  a vector  by  a Toeplitz  matrix 
(or  a partition  of  a Toeplitz  matrix)  has  been  described  in 
Section  2.1.  The  matrix  to  be  multiplied  is  either  the  cross 
covariance  matrix  C or  the  transpose  of  the  measurement 
matrix,  H . In  both  cases,  the  dimension  of  the  matrix  to 
be  multiplied,  which  will  be  denoted  here  by  T,  is  m^n  where 
m is  the  dimension  of  the  vector  to  be  estimated  and  n is  the 
dimension  of  the  data  vector.  The  computational  cost  of 
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carrj'ing  out  this  vector-matrix  multiplication  by  standard 
means  is,  of  course,  ran  (scalar)  multiplications  and  m(n-l) 
additions.  When  m is  small  compared  to  n,  an  ordinary  matrix 
multiplication  procedure  can  be  used  without  affecting  the 
overall  efficiency  of  the  estimation  algorithm.  When  m is 
an  appreciable  fraction  of  n,  the  cost  of  this  simple  opera- 
tion may  dominate  the  cost  of  the  data  processing  procedure 
unless  a special-purpose  algorithm  is  used.  Such  an  algorithm 
is  presented  below. 

By  assumption,  the  matrix  to  be  multiplied,  T,  is  a 
Toeplitz  matri.x  or,  when  m<n,  a partition  of  a Toeplitz  m.atrix. 
It  follows  that  the  elements  of  T are  of  the  form 


0 i j < m-1,  0 < k < n-1  (2.2-6) 


where  t,  is 

X, 

the  appropriate  function. 

In 

the 

case  that  m<n. 

the  m.atrix  T 

must  be  extended  so  that 

it  i 

s a 

(square)  Toeplitz 

matrix.  Thi 

s is  accomplished  by  defin 

ing 

the 

nxn  matrix  T by 

rmi  = 

jr^ 

tj._ , ; 0 £ j . £ n-1 

(2.2-9) 

The  multiplication  of  an  arbitrary  vector,  say  b,  by  the 
Toeplitz  matrix  T may  now  be  carried  out,  with  the  first  m 
elements  of  the  result  containing  the  desired  product,  Tb. 
The  problem  then  is  simply  one  of  developing  an  efficient 
procedure  for  calculating  Tb. 

The  means  for  fast  multiplication  of  T follows  from 
the  circulant  theorem  of  Section  2.2.1.  The  matrix  T can  be 
i.mbedded  in  a 2nx2n  circulant  matrix  C,  as  follows.  The 
function  t^  is  extended  to  the  range  -(2n-l)  £ 1 £ (2n-l)  by 
the  definition 
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C, 


I' 


0. 


n-1 

£ = n (2.2- 

n+1  < £ < 2n-l 


and  c p 

— K 


Cry  - for  1 < £ < 2n-l. 
2n-x.  — — 


The  matrix  C with  elenect 


'k-j 


(2.2- 


is  now  a circulant.  .Moreover,  the  circulant  C has  the 
additional  property  that 


(2.2- 


where  b is  an  arbitrary  n-vector  and  £ is  the  zero  vector  cf 
length  n.  The  vector  Tb  is  the  desired  product  vector,  whil 
the  n-vector  d is  discarded.  The  result  of  Section  2.2.1  is 
applied  to  carry  out  the  multiplication  using  fast  rouner 
transforms . 


Define  B and 


to  be  the  transforms 


B = 


(2.2 


ic 


= Fc 


(2.2- 


with  the  elements 


V t 

^ " (^0'  • • • ’ ®2n-l) 

^c  " ("O’  "l’  '^Bn-l) 


(2.2 


-14) 
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len  frora  Eqs.  2.2-o  to  v.  .2-1  the  vector 


C B = >■  n 


0 0 


^3,  . . 

X JL 


''2n-l®2n 


,-x) 


(2.2-15) 


s the  DFT  of  the  product  vector  in  Eq.  2.2-12.  .^n  application 

f the  inverse  Fourier  transform,  F , produces  the  desired 
esult.  The  overall  algorithm  for  fast  matri.x  multiplication 
s diagrammed  in  Fig.  2.2-1.  The  cost  of  the  entire  procedure 
s dominated  by  the  cost  of  three  applications  of  the  Fast 
ourier  Transform  algorithm,  and  is  proportional  to  n log  n. 


Thus,  the  multiplication  of  a vector  by  a Toeplitz 
matrix  is  equivalent  to  the  convolution  of  two  discrete  func- 
tions of  finite  length.  The  algorithm  presented  here  for 
Toeplitz  matrix  multiplication  is,  in  fact,  equivalent  to  a 
well-known  fast  convolution  algorithm,  which  has  been  described 
in  Refs.  9 and  10. 


The 


The  window 

function  is  related 

to  the  way  in 

which  da 

“Z  3. 

i s 

processed , 

and  is  described  in 

the  next  sect 

ion.  For 

3 

general 

discussion 

of  window  functions 

in  frequency 

domain  da 

iX 

r r 0— 

cessing,  see  Ref.  17,  Chapter  3.  Note  that  one  choici 
is  .iimolv 


W = I 


(2.2-15 


where  I is  the  n^n  unit  matrix.  Thus,  the  class  of  windowed 
Toeulitz  matrices  includes  all  Toenlitz  matrices. 


The  DFT  of  the  windowed  Toeplitz  matrix  in  Eq . 2.2-16 


IS 


T = ECr  = EVrT’Xr 


( 2 . 2-19 ) 


2S 


The  work  of  deriving  a useful  elemen t-bj'-element  expression  for 
this  o'.atrix  is  divided  into  two  parts:  first,  obtaining  an 
expression  for  T in  terms  of  transforms  of  diagonal  and  circu- 
lant  matrices  and,  second,  applying  the  circulant  theorem  of 
Section  2.2.1  to  deduce  an  efficient  computational  formula. 

The  details  of  this  derivation,  which  is  central  to  the  present 
technique,  are  provided  in  Appendix  A. 


The  result  is  the  expression 

2n-l 


J. 

F-  1 = o - 

1=0 


(2.2-20) 


where  the  2n  dimensional  vectors  T and  t.  2.re  defined  as 
follows.  The  complex  vector  H is  the  2n  dimensional  transform 


w 

(20,  1 ’ ■ ' 

^ .2^ 

^2n 

0 

(2.2-21: 


of  the  extended  window  function  where 


H - ("O’  "i "n-l 


( 2 . 2-22  ) 


The  real  vector  _t  is  the  power  spectruri  corresponding  to  the 
extended  circulant  function  c^^  defined  in  Eq . 2.2-10.  That 
is,  2 transform 

T 


- ( = F c 

•0’  -1’  2n-l/  2n  - 


(2.2-23) 


* The  notation  2^  represents  the  complex  conjugate  of  the 

auantitv  T , . 

Jo 

**The  2n  dim.ensional  DFT  matrix,  ^2^^  • is  defined  by  Eq . 2.2-3 
with  n replaced  by  2n . 
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Equation  2.2-20  displays  the  elements  of  as 
weighted  sums  of  a power  spectrum  t_,  where  the  weights  are 
determined  by  a spectral  window  n.  Since  T can  be  internreted 

“ Uj 

as  a frequency  domain  covariance  matrix,  Eq . 2.2.20  is  analo- 
gous to  the  integral  formulas  which  hold  in  the  spectral  repre- 
sentation of  a stationary  stochastic  process  (see,  for  e.xample, 
Ref.  11,  p.  200).  The  application  of  this  formula  for  to 
the  solution  of  the  systems  of  equations  encountered  in  gravity 
quantity  estimation  is  described  in  the  next  section. 

2.3  EFFICIENT  SOLUTION  OF  COVARIAxNCE  EQUATIONS 

The  requirement  in  gravitj'  quantity  estimation  for 
an  algorithm  to  solve  large  systems  of  equations  with  coeffi- 
cients given  by  covariance  m.atrices  has  been  discussed  in 
Section  2.1.  The  solution  of  such  equation  sets  is  tne  single 
most  expensive  operation  in  the  estim.ation  process.  For 
proble.m.s  involving  large  a*mounts  of  data,  the  use  of  standard 
algorith.ms  for  solving  the  covariance  equations  results  in 
prohibitive  computer  time  requirements.  For  gravity  esti.ma- 
tion  problems,  or  any  estimation  problemi  involving  stationary 
models,  the  covariance  imatrix  defining  the  syste.m.  of  equations 
has  a highly  redundant,  Toeplitz  structure.  The  exploitation 
of  this  special  mathematical  structure  makes  it  possible  to 
greatly  reduce  the  co.mputer  processing  time  required  to  solve 
the  covariance  equations.  This  increase  in  efficiency  is 
reflected  in  a decrease,  by  orders  of  miagnitude,  in  the  pro- 
cessing time  required  for  the  solution  of  the  overall  esti.ma- 

new  highly  efficient 


r 


algorithm  for  the  solution  of  systems  of  equations  with 
Toeplitz  structure. 

2.3.1  Algorithm  Overview 

The  new  GEOFAST  algorithm  is  based  on  the  transform.a- 
tion  of  the  system  of  covariance  equations  into  the  frequency 
domain,  and  on  the  use  of  an  appropriate  data  window  to  control 
the  frequency  domain  matrix  structure.  Consider  the  solution 
of  a svstemi  of  eauations 


Tu  = z 


V ^ ; 


where  T is  an  n>^n  Toeplitz  matrix.  For  the  gravity  quantity 
estimation  problem.,  T represents  the  modeled  data  covariance 
miatrix  and  z is  the  data  vector.  The  vector  u represents  the 
solution  vector  required  for  subsequent,  stages  of  processing. 
Note  that  Eq . 2.3-1  m.ay  be  solved  using  any  convenient  coordi- 
nate transf  orm.at  ions . If  a nonsingular  transformation.  A,  is 
appropriately  applied  to  u and  z the  equivalent  system. 


u “ Z 


(2 . 3-2a) 


IS  cctainea.  wnere 


(A  )"^  u 


(2.3-2b) 


Once  Eq.  2.3-2  is  solved,  the  solution  vector  u is  sim.ply 
found  from. 


u = A u 


(2.3-3) 
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Of  course,  the  objective  in  carrying  out  such  a transformation 
is  to  obtain  a system  of  equations.  Eq.  2.3-2,  which  is 
inexpensive  to  solve. 


GEOFAST  employs  a transformation  which  is  composed, 
first,  of  multiplication  by  a data  window  and,  second,  of 
application  of  the  DFT . Multiplication  by  a data  window  is 
represented  by  an  n^n  diagonal  matrix,  defined  by  Ec.  2.2-17. 
The  transf ormation  A under  consideration  here  is  thus 


A = FV  (2.3-4) 

where -F  is  the  DFT  matrix  of  Eq . 2.2-3.  vrith  this  choice  T' 
is  identical  to  the  matrix  T defined  in  Ec . 2.2-19.  The 

Ll; 

m.ain  advantage  of  this  transformation  is  that  it  can  be  made 
to  give  the  matrix  T'  a very  sparse  structure,  with  only 
negligibly  small  elements  outside  a narrow  band  about  the 
main  diagonal.  This  special  structure  is  discussed  in  detail 
in  Section  2.3.2.  A second  imipcrtant  advantage  of  the 
\vindowed  DFT  is  that  it  is  very  inexpensive  to  apply;  window- 
ing itself  is  a very  simple  operation,  while  the  DFT  can  be 
carried  out  with  the  Fast  Fourier  Transform  algorithm. 

A block  diagram  for  the  overall  algorithm  is  presented 
in  Fig.  2.3-1.  The  transformation  described  by  Eqs.  2.3-3  and 
2.3-4  is  used  to  produce  a system;  of  equations  for  which  the 


matrix  T’  is  very 

near 

ly  band-diagonal 

The  m.a 

trix  T 

' i s 

then  approximiated 

by  a 

m.atrix  T.  which 

is  exact 

ly  band-diago 

and  which  differs 

from 

T'  only  in  that 

it  sets 

sm.all  out -of- 

band  elements  to 

zero  . 

The  system,  of 

ecuat ions 

, Ec  . 

2 . 3-2  , 

is  then  solved  with  T'  replaced  by  T^.  . The  matrix  T.  has  the 
band-diagonal  structure 
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where  all'  the  elements  of  T.  are  zero  excent  for  those  in  a 

0 

diagonal  band  consisting  of  the  main  diagonal  and  sub-  and 

superdiagonals.  The  difference  between  T'  and  T , consisting 

c 

of  the  out-of-band  elements  of  T' , can  be  made  as  small  as 
desire^  b3-  adjusting  the  data  window  and  selecting  an  appro- 
priate .matrix  bandwidth,  The  relationship  between  the 

data  window  and  the  size  of  the  neglected  elements  is  explored 
in  Section  2.3.2. 


SOL^JTlO^ 

ij 


Figure  2.3-1 


Overview  of  GiOFAST  Algorithm 


By  applying  the  windowed  DFT , the  system  of  equations 

to  be  solved  has  been  reduced  to  one  with  the  band-diagonal 

matrix  T.  of  Eq . 2.3-5.  This  is  extremely  advantageous  because 

a band-diagonal  sj’stem  of  this  type,  where  .M,.,  is  much  smaller 

than  the  dimension,  n,  of  T - , can  be  solved  in  far  fewer 

numerical  operations  than  a full  system.  Application  of  the 

band-diagonal  implementation  of  Cholesky  decomposition  (Ref.  12) 

2 

requires  numerical  ooerations  prooortional  to  only  M^^n , as 
3 ' ' B 

compared  to  n for  a full  system.  The  results  presented  in 

Section  2 . 3 .'4  indicate  that  values  of  Mg  on  the  order  of  10  miay 

be  adequate  for  T » . For  M_  = 10,  the  resulting  reduction  in 

computer  processing  time  is  proportional  to  (n/lO)*^,  which 

amount’s  to  several  orders  of  magnitude  for  large  gravity 

quantity  estimation  problems.  A consequence  of  the  reduction 

of  the  system  of  equations,  Eq.  2.3-1,  to  band-diagonal  form, 

is  thus  a substantial  reduction  in  the  amount  of  computer  time 

required  to  obtain  a solution. 

The  substitution  of  the  band-diagonal  matrix  T-  for 
the  exactly  transformed  matrix  T'  involves  an  approxim.at on 
which,  of  course,  effects  the  solution  vector  u.  The  m.agni- 
tude  of  this  error  is  controlled  by  the  choice  of  data  window, 
the  number  of  diagonal  bands  included  in  T. , and  the  size  of 
a damping  factor  introduced  in  the  Cholesky  decomposition 
process  (Section  2.3.3).  By  appropriate  choice  of  these  con- 
trolling parameters,  the  approxim.at ion  error  in  u can  be  made 
as  small  as  desired.  It  is  to  be  emphasized  that  what  is  being 
considered  here  is  computational  error  in  approximating  T'  with 
T . , rather  than  statistical  error  in  computing  an  estim.ate  from 
datai  The  statistical  error  is,  of  course,  subject  to  the  usual 
considerations,  and  not  controlled  in  any  way  by  the  computa- 
tional parameters  discussed  here.  The  com.putational  error  in 
the  new  algorithm,  and  mechanism.s  for  controlling  it  are 
analyzed  in  Section  2.3.3. 
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2.3.2  Computation  of  the  Transformed  Covariance  Matri.x 


The  computation  of  the  transformed  Toeplitz  matrix, 

T. , is  a major  portion  of  the  new  algorithm.  This  part  of  the 
algorithm  is  particularly  important  for  two  reasons: 

• The  expression  for  T5  shows  how  the  data 
window  determines  matrix  bandwidth. 

• Special  computational  methods  are  required 
to  keep  the  computer  time  required  for 
generating  the  ele.ments  of  T.  to  a minimum. 

^ C 

Both  reasons  m.ake  the  form,ulas  for  the  elements  of  T.  and  the 
manner^  in  which  they  are  imiplemented  central  to  the  overall 
algorithm . 


As  described  in  Section  2.3.1,  the  matrix  T.  is 

0 

defined  conceptually  by  applying  a linear  transf orm.at ion  given 
in  Eq . 2.3-4  to  the  matrix  T to  obtain  the  matrix  T'  from. 

Eq . 2.3-3,  and  setting  to  zero  all  elements  of  T'  outside  a 
predeterm.ined  bandwidth.  The  word  "conceptually”  is  used 
here  because  the  detailed  algorithm,  explained  below,  modifies 
the  DFT  .matrix  E shown  in  Eq . 2.3-4  to  use  onl3'  half  of  the 
frequencies  present.  It  will  become  apparent  that  this  modi- 
fication simply  eliminates  the  redundancy  which  occurs  because 
the  data  are  always  real-valued.  An  expression  for  T',  before 
any  modification  is  made  to  F,  is  available  directly  from 
Section  2.2.3.  Since  T ’ =T  Eq . 2.2-20  provides 


2r-l 


where  the  subscripts  of  the  elements  of  T'  have  been  arranged 
so  that  the  elements 

T.l  ^ £ n-m-1  (2.3-7) 

- ' -Ik  , k+m 

T h 

represent  the  m ' diagonal  band.  The  functions  and  are 
defined  by  Eq . 2.2-21  and  2.2-23  as  described  in  Section  2.2.3. 

The  impact  of  the  window  function  on  the  structure 

- K 

of  T'  can  be  seen  from  Eg.  2.3-6.  The  elements  of  the  m'"” 
diagonal  band  are  obtained  by  convolving  the  power  spectrum 
T,  with  the  weighting  function. 


(It)  = 2, 


£-2m 


(2.3-8) 


which  consists  of  the  transformed  data  window  multiplied  by 
its  conjugate  shifted  2m,  elements  to  the  right. 


Two  typical  window  functions  and  their  transform.s 
are  shown  in  Fig.  2.3-2.  These  windows  are  special  cases  of 
the  Kaiser  window  function  which  has  proved  very  useful  in 
signal  processing  applications,  and  is  defined  in  Refs.  17 
and  18.  Transforms  of  data  windows  are  characterized  by  a 
m.ainlobe  and  a series  of  sidelobes  of  considerably  reduced 
levels.  As  can  be  seen  from  Fig.  2.3-2,  the  disparity  between 
mainlobe  and  sidelobe  levels  becomes  more  pronounced  as  the 
width  of  the  m.ainlobe  is  increased.  This  behavior  is  reflected 
in  the  num.erical  structure  of  the  matrix  T'.  The  elements  of 
T'  corresponding  to  the  shaded  area  of  Fig.  2.3-3  cannot  be 
neglected  because  the  mainlobes  of  2,  and  2.  r,  overlap  for 
these  elements.  However,  the  elements  of  T'  outside  the 
shaded  area  of  the  figure  a'-e  the  product  of  the  sidelobes 
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of  r. , together  with  the  power  spectrum  t,.  When  the  sidelobes 

X Aw 

are  made  sufficiently'  small,  these  elements  may  be  neglected. 
The  key  idea  in  the  band-diagonal  algorithm  is  the  control  of 
out-of-band  covariance  elements  by  choice  of  a data  window 
with  sufficiently  small  sidelobes. 

R-30299 
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Figure  2.3-3  Structure  of  Covariance  Matrix 

Under  Fourier  Transformation 

A final  transformation  must  be  applied  to  form  the 
frequency  domain  covariance  matrix  into  a band-diagonal  struc- 
t'ure.  This  transformation  corresponds  to  deleting  from  the 
data  redundant  components  at  negative  frequencies.  The  nega- 
tive freauencv  com.Donents  are  redundant  in  a statistical 


sense  because  it  is  known  that  the  data  niust  always  be  real- 
x-alued.  Put  another  way,  the  sine-cosine  transform  coeffi- 
cients are  sufficient  to  describe  the  data,  rather  than  the 
full,  complex  set  Of  Fourier  coefficients.  The  practical 


importance  of  this  last  transformation  is  the  elimination  of 
the  elements  in  the  upper  right  and  lower  left  corners  of  T’  as 
shown  in  Fig.  2.3-3.  With  the  elimination  of  these  elements,  a 


The  firsi  n/2-i-l  rows  of  H generate  the  cosine  coefficients 
of  the  data,  while  the  last  n/2-1  rows  generate  the  sine 
coefficients.  The  result  of  applj'ing  the  modified  transfor- 
.mation,  Eq . 2.3-9,  to  the  definition  of  T'  is  shown  in  Fig. 
2.3-4.  The  matrix  T'  is  now  a simple  band-diagonal  matrix, 
neglecting  the  small  elements  due  to  the  sidelobes  of  f.,.  The 
sine-cosine  transformation  has  two  additional  advantages: 


• The  nonzero  elements  of  T'  may  be  par- 
titioned into  blocks  of  dimension  n/2*l 
and  n/2-1  corresponding  to  cosine  and  sine 

. coefficients.  The  dimension  of  the  largest 
system  of  equations  to  be  solved  is  thus 
approximately  halved. 

• The  elements  of  I'  are  real-valued 
eliminating  the  need  for  complex-valued 
numerical  operations. 


Other  than  these  modifications,  the  banded  structure  of  the 
transf or.Tied  covariance  matrix  is  unaffected  by  the  introductio 
of  H . 


Aleorithm.  Perform.ance 


The  use  of  an  appropriate  data  window  gives  the  DF7 
of  the  covariance  matrix  a numerical  structure  which  m.ay  be 
readily  exploited  in  solving  equations.  The  analysis  of 
Section  2.3.2  shows  that  T.  will  have  the  band-diaaonal  form 
shown  in  Eq . 2.3-5  and  Fig.  2.3-4.  The  m.atrix  T.  differs 
from,  the  exact  frequency  dom.ain  covariance  m.atrix  T’  only 
in  that  it  neglects  the  very  sm.all  elements  which  fall  out- 
side the  diagonal  band.  A very  efficient  approxim.ate  solu- 
tion to  the  covariance  equations  is  found  by  replacing  the 
exact  m.atrix  T'  by  T.  in  Eq.  2.3-2.  The  special  structure  cf 
T is  exploited  by  employing  the  band-diagonal  implementation 
of  the  Cholesky  decomposition  algorithm.  (Ref.  12.  Section  2.3) 
for  the  solution  of  linear  equations.  The  resulting  error  in 


40 


corr.'3_emen“ary  cnan&e  in  xns  naia  vi'inaow  runction  as  csscnoe 
belo'-v.  Bounds  cn  the  con'.tutational  error  arising  fror.  this 
source  are  presented  here. 

The  computational  details  of  the  handed  Cholesky 
decomposition  algorithm  are  provided  in  Ref.  12.  Only  two 
features  of  the  algorithm  are  important  to  the  present 
application : 

• The  banded  Cholesky  algorithm  is  an  exact 

procedure  (within  numerical  error)  for  solving 
a system,  of  equations  with  banded  structure 
of  the  form,  of  Eq . 2.3-5. 


• The  number  of  numerical  operations  required 

o' 

is  approximately  C(Mg+l)“n,  where  C is  a 

sm.all  constant  dependent  on  the  particular 
im.plementation . 


For  sm.all  values  of  .VL 

D 

solving  Eq.  2.3-2  with 


(sav  .\I_  < 25)  the  numerical  work 
the  banded  Cholesky  procedure  is 


of 

actually 


less  than  that  required  to  compute  the  elements  of  T.. 

0 


The  error  bounds  given  in  this  section  are  parameter- 
ized by  a dam.ping  factor,  c,  and  a resulting  quantity  which 
will  be  called  data  de-em.phasis . In  order  to  guarantee  that 
the  error  in  replacing  T'  by  T.  in  Eq . 2.3-2  produces  a sm.all 
error  in  the  solution,  it  is  necessarv  to  m.odifv  T'  bv  intro- 

f ... 

ducing  a s.m.all  amount  of  ficticious,  white  noise.  The  system, 
of  equations,  Eq.  2.3-2,  is  replaced  by 


(T-  e cl)u'  = z'  (2.3-11) 

When  Eq.  2.3-11  is  solved  by  Cholesky  decom.pcsition  T'  is 

aoDroxim.ated  bv  T.  and  Ea . 2.3-11  becomes 
* * * 0 

(T^  - cDu"  = z'  (2.3-12) 


The  num.ber  d is  a sm.all  constant  which  can  be  chosen  to  trade 


off  co.T.putat ional  effort  for  com.putat ional  error.  The  modifi- 
cation reoresented  bv  Ecs.  2.3-11  and  2.3-12  is  necessarv 


oecause  c: 


eigenvalues  introduced  in  T’  bv  the 


use  of  highlv  taoered  data  windows.  It  can  be  shown  That  the 


square  of  the  sm.allest  weight  in  the  data  window  is  an  upper 

bound  to  the  sm.allest  eigenvalue  of  T'.  These  sm.all  eigen- 

_ 1 

values  introduce  very  large  eigenvalues  in  (T')  *,  which 


greatly  am.plify  the  error  (u'-u")  due  to  the  difference 
between  T'  and  7-.  That  is,  the  system  of  Eq . 2.3-2  is 
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ill-conditioned*  when  windows  are  used.  By  contrast,  the 
matrices  T’  +61  and  + cl  can  have  no  eigenvalues  smaller 

w 

than  6 so  that  Eqs . 2.3-11  and  2.3-12  are  well-conditioned. 

The  actual  impact  of  this  modification  on  the  error  (u'-u") 
is  reflected  in  the  error  bounds  presented  in  this  section. 

The  modification  of  T'  represented  by  Eq . 2.3-11 
has  an  important  interpretation  in  terms  of  the  original 
estimation  problem..  The  original  system  of  equations  involv- 
ing T',  Eq . 2.3-2,  is  simply  the  transform,  of  the  covariance 
equations.  Ec . 2.3-1.  By  inverting  the  transformation  pro- 
cess, the  syste.m  of  covariance  equations  corresponding  to 
Eq . 2.3-11  is  found  to  be 

(T  D)  u = z (2.3-13) 

where  D is  the  diagonal  matrix 


D = 

5 

diag 

^ -2 
(Wq  , 

-2 

^1  ’ 

-0 

(2 

. 3-14 

Equation 

2 . 

3-13 

shows  that 

the  addition  of 

a multiple 

r\  * 

the 

identity 

m.a 

tr  ix 

to  T ' 

corr 

esponds  to  the  ad 

dition  of 

the 

sam.  e 

.multiple 

of 

the 

inverse  of 

the  window  m.atri 

X squared 

to  T 

This  m,ay 

be 

inte 

rpreted  as 

a model  for  the 

addition  c 

f un 

- 

correlate 

d 

measu 

rement  noi 

se  to  the  data. 

(Note  that 

no 

noise 

is  actually  added  to  the  data. ) The  variance  of  this  addi- 
tional noise  changes  with  loca'^ion  along  the  data  track  in 
proportion  to  the  square  of  the  inverse  of  the  window  function. 
The  variance  of  the  added  noise,  which  is  large  near  the  ends 
of  the  track,  and  small  near  the  center,  is  called  the  data 
de-e.m.phasis  function  and  is  displayed  in  Eig.  2.3-5  for  a 
tvnical  Kaiser  data  window. 


’•'Ill-conditioning  is  discussed  in  Ref.  20  and  Appendi.x  3. 
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normalized  data  length 


NORMALIZED  DATA  LENGTH 


zigure  2.3-5  Data  Window  and  Resulting  De-r.mphasis 

(Kaiser  Window,  ^2=6,  6=10“^,  n=12S) 

The  modification  represented  by  Eqs . 2.3-11  and 
2.3-13  changes  rhe  esti.mation  problem  which  is  actually 
solved  by  the  nu.merical  algori'^hm.  In  the  modified  problem 
the  quality  of  the  data  is  gradually  de-emphas ized  as  the 
ends  of  the  da~a  track  are  approached.  Points  along  xhe 
track,  for  which  the  signal-to-noise  ratio 


S 


k 


(2.3-15) 


is  less  than  one,  may  be  regarded  as  having  been  lost  due  to 
the  introduction  of  computational  data  de-emphasis.  The  frac- 
tional number  of  data  points  for  which  is  less  than  one 
defines  a perform.ance  index,  the  data  de-emiphasis  index,  which 
is  associated  with  a given  dam.ping  factor,  5,  and  a given 
window  function.  The  data  de-emphasis  resulting  from  a number 
of  com.binat ions  of  dam.ping  factor  and  m.atrix  bandwidth  is 
presented  in  Fig.  2.3-6. 


f 


Figure  2.3-6  Data  De-Emphasis  Resulting 

from  Computational  Damping 


An  important  concern  in  the  use  of 
method  to  compute  estimates  from,  data  is  the 
in  esti.mation  error  due  to  the  computational 
case  rf  the  band-diagonal  method  fcr  solving 
tions,  computational  error  in  the  estim.ate  a 
error  vector 


an  annroxim.ate 


expected  increas 
error.  In  the 
covariance  equa- 
rises  because  the 


:u  = u'  - u" 


(2.3-1 


is,  in  general,  nonzero.  Since  the  data  vector  z is  a random, 
quantity,  cu  and  the  resulting  increase  in  esti.mation  error, 
fx,  are  also  rando.m.  The  m.agnitude  of  this  error  can  be 
measured  by  the  ratio 


E(cx‘  ox) 

ECx^  X) 


(2.3-1 


mV  r 


ienotes  the  ensemble  expectation  operator. 


which  is  the  relative  rms  computational  error  in  x.  A realistic 
upper  bound,  w'hich  may  be  readily  com.puted,  for  the  quantity  c 
is  derived  in  Appendix  B.  This  bound,  as  well  as  the  actual 
value  of  c , is  a function  of  the  damping  factor  c and  the 
number  of  superdiagonal  bands.  .M„ , used  in  T-.  Figure  2.3-7 
presents  a set  of  values  for  the  error  bound  obtained  when  T 
is  generated  by  a third  order  Markov  model  for  the  gravity 
anomaly  (Ref.  13). 

Figures  2.3-6  and  2.3-7  show  the  two  major  measures 

of  performance,  data  de-emphasis  and  computational  error,  as  a 

function  of  matrix  bandwidth  , which  in  turn  determines  the 

o 

amount  of  computation  required.  For  example,  winh  seven  suner- 

lu 

diagonal  bands  and  a damping  factor  of  10  * the  error  bound  is 
O.lFc.  The  resulting  data  de-emphasis  is  about  22'5c.  It  is 
possible  to  achieve  the  sa.me  error  bound  with  a damping  factor 
of  10”  by  using  nine  superdiagonal  bands  and  performing  more 
com.putation . In  this  case  the  data  de-em.phasis  is  reduced  to 
17’^c.  There  is  thus  a tradeoff,  controlled  by  the  damping 
factor,  between  data  de-emphasis  and  comiputational  effort.  A 
given  level  of  com.putat ional  accuracy  m.ay  be  obtained  with  as 
sm.all  a level  of  data  de-emphasis  as  desired,  at  the  expense 
of  a larger  m.atrr.x  bandwidth  and  increased  computing  time. 

The  effect  on  com.putat  ional  error  and  data  de-emphasis  of 
reducing  the  damping  factor  is  demonstrated  in  Fig.  2 . 3-S . 

2.4  APPLICATION  DISCUSSION 

To  provide  perspective  on  the  use  of  GFOFAST,  a sample 
application  to  gravity  gradiometer  data  is  outlined.  For  simp- 
licity the  example  treats  only  one  of  the  five  independent 
elements  of  the  gravity  gradient  tensor.  In  most  actual  app- 
lications, however,  all  of  the  measured  gradients  would  be 
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Figure  2.3-7  Helaxive  RMS  Compu'ax ional  Error 

for  Band  Diagonal  Method 

R-30297 


PERCENT  DATA  DE-EMPHASIS 

Figure  2.3-8  CoTiputational  Error  and  Data  De-Emphasis 

Regimes  for  Various  Damping  Factors 
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processed  as  rnultisensor  measurements.  Extension  to  many  sen- 
sors is  straightforward  (Ref.  1). 

Suppose  that  10,000  measurements  of  the  anomalous,  along- 
track,  along-track  gravity  gradient  are  available  at  0.1  km.  inter- 
vals. i.e.  on  a 1000  km  track.  Such  measurements  could  be  the  out 
puts  of  an  airborne  gradiometric  survey  which  have  the  reference 
spheroid  gravity  gradient  field  removed  and  have  been  pre- 
filtered (e.g.  averaged).  Opti.mal  estimates  of  the  vertical 
deflection,  at  the  earth’s  surface  below  the  data  track  are 
desired  at  a spacing  corresponding  to  that  of  the  airborne 
measurements . 


The  optimial  estimiate  of  I is  given  by  Zq . 2.1-2  where 


X = 


(2.4-1) 


z = 


V 


(2.4-2) 


where  v is  noise.  In  Eqs . 2.4-1  and  2.4-2,  g is  gravity,  t de- 
notes the  along-track  direction  and  2--  -s  the  vector  of  anomalous 
along-track,  along-track  self-gradients.  The  covariance  .m.atrices  of 
Eq.  2.1-2  must  be  specified  by  a suitable  statistical  model.  The 
attenuated  white  noise  (A\VN)  model,  described  in  Refs.  1 and  21. 
is  appropriate  for  computing  the  required  covariances.  Using  the 
AWX  model,  the  autocovariance  of  the  along-track  gradient, 

2^2 <^t),  is  given  by 


P.  9 "-2  i 

^ — t 9C-n°-  ^ \ 1 

[T^-8(D+h)^] 

+ 1 

[t^  - 4(D-h)^] 

1 -1-  D 

[-2  _ 4(D+h)“] 

11/2 

(2.4-3) 


where  h is  the  altitude  at  which  the  gradient  measurements  are 
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taken,  t is  the  shift  distance  and  R is  the  measurement  error 
covariance  matrix.  The  model  parameters,  a__  (the  rms  along- 
track,  along-track  gradient)  and  D (a  measure  of  the  rate-of- 
decrease  of  the  autocovariance  function),  are  determined  by  statis- 
tical analysis  of  the  data  as  described  in  Ref.  21.  Typical  val- 

* 

ues  to  be  expected  are  24  EU  and  35  km  respectively.  The  cross 
covariance  between  the  anomalous,  along-track,  along-track  self- 
gradient at  altitude,  h,  and  the  along-track  gravity  disturbance 
on  the  ground  is  given  by  the  AWN  model  as 


r 


xz 


(t)  = 


( 2D-h ) 


1 

t i 1 

[T^-2(2D-h)^l 

1 + ' 

r 2 2- 

L-  -^(2D^h)  ! 

V 

i 

[t^  - (2D  . h)"] 

11/2 

The  elements  of  the 

jk  " '"zz  ^'jk^ 

[^xz' 

jk  ■ ^xz  ■ jk- 

r 

= ( h - j ) _ 

zz 


xz 


(2.4-4) 


are  given  dv 


(2.4-5) 

(2.4-6) 

. k = 1,  2 ...  10,000  (2.4-7) 


and  1 is  the  data  spacing,  0.1  km.  The  indices  j and  k are 
data  ele.ment  counters. 


In  principle,  Eq.  2.1-2  could  now  be  solved  directly 


:or 


;h9  desired  deflection  estimates.  However,  in  practice, 
the  10,000  X 10,000  dimensionality  of  would  make  the  C 


In  addition. 


zz 

remult i- 


operation  computationaly  prohibitive 

olication  of  C~~  bv  C vould  be  exceedinglv  burdensome, 
zz  • xz  ® ■ 

The  initial  operation  of  the  GEOrAST  algorithm  is  to 
multiply  the  data  vector,  element-by-elemer.  c , by  the  Kaiser 


_o  _o 

*1.0  EU  = 0.1  mgal/km  = 10  ''  sec 


AO 


W 


X 

window  exemplified  in  the  lower  right  graph  of  Fig.  2.3-2. 

Note  that  in  the  windowing  operation,  as  with  all  of  the 
operations  involved  in  GEOFAST,  the  count  of  scalar  arithmetic 
.procedures  is  of  the  order  of  the  number  of  elements  in  the  data 
vector  (except  FFTs  which  go  as  n log  n).  Following  windowing, 
the  gradiometer  data  is  fast  Fourier  transformed  into  the 
"intermediate"  frequency  domain  data  vector,  z. 

The  transform  of  the  Kaiser  window  function,  1,  is 
computed  by  Eq . 2.2-21,  an  operation  involving  a FFT  and  scaling. 
The  "window  shape"  variable,  Mg,  is  chosen  from  Fig.  2.3-7  tc 
provide  acceptable  accuracj-  and  reasonable  computation  time. 

Note  that  the  choice  of  Mg  specifies  the  number  of  bands  whtch 
are  ultimately  retained  in  the  matrix  inverse  computation. 

V.'hen  Mg  is  selected,  the  value  of  the  data  de-emphasis  damping 

factor,  c,  is  also  chosen  from  Fig.  2.3-S.  Transf orm.ation  of. 

n,  by  the  miatrix  H (given  by  Eq . 2.3-10)  follows,  i.e. 

£ g = K ^ (2.4-S) 

The  circulant  vector  c (first  column  of  the  circulant 
m.atrix,  C)  is  formed  using  Eq . 2.2-10.  The  first  10,000  elements. 

c^,  are  given  by  the  AWN  model  along-track  self-gradient  covar- 

iance function  of  Sqs.  2.4-3  and  2.4-5  with  j=l  and  k = 1,  2 ... 
10,000.  The  FFT  of  c is  computed  using  Eq . 2.2-23  which  provides 
the  elements  of  the  gradient  auto-spectrum  t_.  Equation  2.3-6 
is  now  invoked  to  comnute  the  M_  diagonal  band  elements, 

D 

FT-I  , 1.  4.  . which  will  be  supplied  to  the  inverse  calculation. 

I cl  ic , .i  ^ m 

In  the  application  of  Eq.  2.3-5  the  0..  elements  used  are  taken 


’'In  Fig.  2.3-2  only  half  of  the  window  is  depicted.  In  actual 
fact  , the  Kaiser  window  is  sj'metric  over  the  data  interval  and 
weights  data  in  the  center  most  hea\’ily. 
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from  the  transformed  window  function  vector,  To  the  main 

diagonal  elements  of  T . , the  value  of  the  data  de-emphasis 
factor,  0,  chosen  at  the  outset,  is  now  added. 

The  equation  set 

T.  u'  = 2'  (2.4-9) 

C 

where  T-  is  banded  diagonal  and  z'  the  windowed,  transformed 
gradient  data  vector,  is  now  completely  specified.  Application 
of  the  banded-Cholesky  decomposition  routine  to  Eq . (2.4-9) 

r 

provides  the  "intermediate  solution"  vector  u . 

To  prepare  the  vector  u'  for  the  fast  multiplication 
portion  of  the  estimation  process  (i.e.,  accounting  for  the 
elements  of  the  cross  correlation  between  gradients  and  de- 
flections. C ) the  transformation  given  bv  Ecs.  2.3-3  and  2.3-9 
must  be  applied,  namely, 

b = vrF  H u (2.4-10) 


i 


i 


I 


1 

I 


Since  the  window  matrix, 
sparse  structure  sim.ilar 
implemented  without  adve 


V.’ , is  diagonal  and  H has  a 

to  that  of  H,  Eq . 2.4-10  is  easily 

se  effect  on  the  efficiency  of  GEOEAST. 


The  vector,  b,  is  the  solution  to  the  inverse  portion 


0 ^ 

the 

proble.m 

E ast 

multiplication 

is  now  used  to 

obtain  the 

mai 

ri.x 

product 

^.xz 

The  gradient- 

deflection  cross 

covariance 

r i V 

deficit 

ions  of 

Eqs.  2.4-3  and 

2.4-6  are  used 

provide  elv 

ments  for  the  circulant  vector,  c.  of  Eq . 2.2-10.  Equation  2.2-10 
is  applied  to  in  the  same  fashion  as  previously  accomplished 

for  in  the  inversion  portion  of  GE0F.A,ST.  The  vector  c 

undergoes  the  EFT  operation  as  indicated  in  Eq . 2.2-13a  to 
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provide  The  solution  vector,  b,  to  the  inverse  portion 

of  the  problem  is  augmented  with  zeros  and  fast  Fourier  trans- 
formed b3’  Eq.  2.2-13b.  The  result,  B,  is  multiplied  element- 
by-element  by  and  scaled  as  indicated  in  Eq.  (2.2-15).  This 
result  is  inverse  Fourier  transformed  to  complete  the  esti.ma- 
tion  procedure.  The  first  10,000  elements  of  the  inverse  trans- 
form are  the  desired  along-track  vertical  deflection  estimates, 

/s 

Z-  It  is  appropriate  to  emphasize  that  the  deflection  estima- 
tion problem  outlined  above  was  chosen  only  for  conceptual 
clarity.  The  actual  solution  of  this  problem  would  not  pose  a 
severe  test  of  GEOFAST’ s capabilities. 
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3.  CONCLUSIONS  AND  RECOMMENDATIONS 

y In  Ref.  1,  TASC  derived  a frequency  domain  algorithm 

for  efficiently  processing  gravimetric  sensor  data.  This 
report  presents  significant  alterations  to  the  algorith.m  based 
upon  e.xplicit  compensation  for  the  adverse  computational  effects 
resulting  from  the  finite  length  of  gravimetric  data  records. 
Using  the  new  algorithm  the  accuracy  of  gravity  quantity  esti- 
mates is  improved  without  noticable  sacrifices  in  computer 
speed.  As  a result,  the  enhanced  algorith.m,  GEOFAST,  is 
applicable  to  a very  broad  class  of  gravity  quantity  estima- 
tion problems. 

The  keynote  of  GEOFAST  is  judicious  application  of 
fast  Fourier  transform  techniques.  However,  since  these  tech- 
niques exploit  certain  m.athematical  structures,  the  gravity  data 
and  the  associated  covariance  models  must  exhibit  these 
structural  forms.  The  requirements  are 

• Input  data  must  be  regularly  spaced  on  ? 

Cartesian  grid.  If  the  original  data 
not  rectangularly  gridded  it  can  be  con 
verted  by  an  appropriate  preprocessing 
stage  involving  averaging  or  interpolation. 

• Stationary,  self-consistent  statistical 
gravimetric  models  must  be  used. 

Currently  the  algorithm  is  formulated  in  terms  of 
estimation  along  and  above  a one-dimensional  survey  track. 
Preliminary  investigation  suggests  that  extension  to  a two- 
dimensional  region  can  be  accomplished  without  fundamental 
difficulty  but  will  involve  considerable  attention  to  mathe- 
matical details.  Prior  to  such  extension,  it  is  recommended 
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that  a range  of  tests  be  conducted  which  exercise  the  algo- 
rithm with  both  simulated  and  real  gravimetric  data.  The 
purpose  of  the  tests  is  to 

• Validate  the  algorithm  by  performing 
several  "end-to-end"  checks 

• Establish  performance  for  very  large 
data  samples 

• Verify  theoretical  estimates  of  accuracy 
and  com,puti*r  speed 

• Identify  possible  sources  of  numerical 
ill-conditioning  or  dynam.ic  range 
lim.it  at  ion . 

The  new  algorithm.,  when  fully  tested,  can  be  expected  to  pro- 
vide a m.ature  tool  for  esti.m.ating  gravity  field  quantities 
from.  data. 
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APPENDIX  A 


TEE  TRANSFORM  OF  A yiNDOyED  TOZPLITZ  MATRIX 

Consider  an  n^n  Toeplitz  matrix  T with  elements 

[T]jj.  = 0 < j,  k < n-1  (A-1) 

corresponding  to  the  function  t-  where  -(n-1)  5 ^ ^ (n-1). 
Let  the  window  function  w.  be  given  for  C < 1 < n-1,  and 

X — — 

define,  the  n^n  matrix 


w = diag  (Wq,  w^ 


(A-2) 


A windowed.  Toeplitz  matrix  is  a matrix  of  the  form. 
V.'T’A' . The  Discrete  Fourier  Transform.  (DFT)  of  such  a m.atrix 
is  given  bp 


T = F TTWF 
X n n 


(A-3) 


wnere 


iy,]  = i exD 


0 < j , k < n-1 


(A-4) 


is  the  n-dimensional  DFT  matrix.  The  purpose  of  this 
appendix  is  to  derive  a computationally  efficient  form^ula 
for  T by  expression  it  in  terms  of  a circulant  matrix. 


*The  symbol  i stands  for  v'-l  and 
transpose . 


indicates  com.plex  conjugate 


A circulant  matrix  (Ref.  8)  is  a Toeplitz  matrix  C 

ith  t , = t The  fundamental  propertv  of  circulant 

— A.  n — X, 

atrices  is  that  they  are  diagonalized  by  the  DFT.  That  is, 

(A-5) 


F CF ' = diag  (d) 
n n — ^ 


d = /n  F 


n — 


(A-6) 


If  C is  a real  symmetric  matrix  then  its  eigenvalues,  dj^^,  are 
real  and  Fas.  A-5  and  A-6  hold  with  F and  F’  interchanged. 
Also  if  d defines  an  arbitrary  diagonal  matrix  then  its  trans- 
form is  a circulant  C.  These  nrouerties  are  used  below. 


The  matrix  T is  extended  to  a 2nx2n  circulant  m^atrix 
T by  the  definitions 


[T]  = c, 

Jk  k-j 

0 <_  j , k _<  2n-l 

(A-7) 

where 

1 'c 

1 ^ 

0 < A £ n-1 

c « = 

X 

I ° 

X = n 

(A-6) 

1 ^i-2n 

< X < 2n— 1 

and  c 

s=  ^ ^ n *»•  1 < 

X 2n-x  — 

i ^ 2n-l. 

The  matrix  T has 

the 

partit 

ioned  formi 

“T  T'  "I 

T = 

i 

T ' T 1 

L ^ _i 

(A-9) 

where 

T'  is  a Toeplitz  matrix  cor 

responding  to  the 

function 

( "i-n 

A > 0 

*C  I “ 

! ° 

A = 0 

(A-10) 

( ^A+n 

A < 0 
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A 2n^2n  extended  version  of  W is  defined  by 


(A- 11) 


The  windowed  Toeplitz  matrix  '\’TW  may  now  be  expressed  in  terms 
of  the  diagonal  matrix  W and  the  circulant  matrix  T by 


nnw 

V 0 


V/TV 


(A-12) 


Introducing  the  nx2n  sampling  matrix  S defined  by 


/l  0 . 

0 0 1 


s = 


(A-13) 


may  be  verified  thal 


T = WTIV  F„„S' 

u)  2n  2n 


(A-14) 


where  F„  is  the  2nx2n  Discrete  Fourier  Transformation  defined 
2n 

by  Eq . A-4  with  n replaced  by  2n.  Equation  A-14  expresses  T 
as  the  result  of  sampling  every  other  row  and  column  of  a 
matrix  which  is  the  transform  of  a product  of  diagonal  window 
matrices  and  a circulant  matrix. 


The  application  of  the  circulant  properties  stated 
above  produces  the  desired  expression  for  T^ . Equation  A-14 
may  be  written  as 
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T = S W T W (A-15) 

ijj  U)  U LJ 

u’here  W and  T are  the  transformed  matrices 

I' 


w 

to 

= ^2n 

^2n 

(A-16) 

T 

U) 

= ^2n 

(A-17) 

Taking  the  coningate  of  Eqs . A-5  and  A-6  shows  that 
the  first  of  these  matrices  is  a circnlant  , wi-ch  the  elements 


I v; 


'k-j 


, 0 < j , k < 2n-l 


(A-IS) 


where  2,.  is  sriven  by  the  transform 

\T 


= o 


X /W  \ 

X 

v'2^ 

the 

window  function 

(A-19) 


'"n-1^ 


( A-20 ) 

It  also  follows  directlv  from  Ecs . A-5  and  A-6  that 


;he  matrix  T of  Eq.  A-17  is  diagonal,  of  the  form 

U) 


- diag  (Tq,  , '2n-l^ 


(A-21) 


where  t = (t,-,  t.  , . . . , ,)  is  defined  bv  the  vector  c of 

— U-L  ^n— 1 ~ 

Eq . A-S  via 


' c 


(A-22) 
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forms 

basic 


The  expression  of  Eq . A-15,  together  with  the  specia 
for  W and  T described  by  Eqs . A-18  and  A-21 , is  the 

U U) 

result  of  this  appendix  in  matrix  form. 


A useful  element-by-element  formula  for  T is  a con- 

sequence  of  Eq . A-15.  Substitution  of  the  definitions  from 

Eqs.  A-13,  A-18,  and  A-21,  with  some  simplification,  produces 
* 

the  expression 


2n-l 


£=0 


"£-2k 


"£ 


(A-23  ) 


where  :r,,  is  the  transform  of  the  window  function,  Eq . A-19, 
and  is  the  power  spectrum  of  T,  defined  by  Eq . A-22 . 
Equation  A-23  provides  a computationally  efficient  formula 
for  the  calculation  of  T . 

OJ 


*The  symbol  .1,  represents  the  complex  conjugate  of  the 
quantity 
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APPENDIX  B 


r 

ESTIMATION  ERROR  BOUNDS 


In  xhis  appendix  an  error  bound  is  derived  for  least 
squares  estimation  using  an  approximate  measurement  covariance 
matrix.  This  bound  is  particularly  useful  in  the  design  and 
analysis  of  fast  estimation  techniques.  Furthermore  the  bound 
is  relati^'ely  easy  to  compute  in  terms  of  known  matrices.  The 
result  is  closel:.'  related  to  well-known  error  bounds  fcr 
deterministic  linear  systems  (Ref.  12,  Section  8 and  Ref.  20, 
Chapter  2).  The  derivation  assumes  that  all  matrices  are  real 
but  readily  extends  to  complex  miatrices  if  "transpose"  is 
replaced  by  "conjugate  transpose." 

Consider  the  least  squares  esti.mate 


where  z is  the  measurement  vector  of  dimension  n,  C is  the 
— zz 

measurement  covariance  matrix,  x is  the  estimate  vector  of 

dimension  m,  and  C is  the  cross-covariance  between  x and  z. 

xz  — — 

An  approximate  covariance  m,atrix 

C.  = C ^ oC  (B-2) 

c zz 

is  constructed  and  used  to  obtain  an  approximate  esti.mate 


It  follows  that 


6x  = -C  (I-C. 

— XZ  0 


'C  )C  z 
zz  zz- 


(3-4) 


where  the  matrix 


R = I 


zz 


(B-5) 


is  called  the  residual  matrix.  Since  x and  z are  random 
variables  an  appropriate  error  measure  is  the  relative  rms 
error  given  by 


E [ 5x ' ox] 
E[x‘x] 


(3-6) 


where  E denotes  the  ensemble  expectation  operator. 


Introducing  the  matrix  trace  operator,  Tr , and  using 


E [zz^]  = C 


zz 


E [XX]  = Tr  E[xx  ] 


the  following  expressions  are  obtained. 

E [e;y£l  = Tr.|c^^Rr^RTcT^} 


E tiT£]  = 


(3-7) 

(3-S) 


(3-9) 


(B-10 


Since  C _ is  positive  definite  it  has  a unique  symmetric 
zz  ^ 

square  root  C~  and  Eq . B-9  may  be  rewritten 


E [iS'ix]  = Tr  |c  C"‘R^C"^C^  } 
— — ( xz  zz  zz  xz } 


(B-11) 


where  the  matrix 


~ +* 
R = C ^RC  “ 
zz  zz 


I - C-  C~^C^ 
zz  zz 


is  a symmetric  form  of  the  residual  matrix, 


(B-i2) 


complete 
Tr  {AB} 


Two  properties  of  the  trace  operator  are  needed  to 
the  derivation.  The  first  is  the  identify 
Tr  {BA;  and  the  second  the  inequality 


Tr  {AB}  < IjAil^  Tr  {B} 


;re  the  euclidean  m.atrix  norm  is  defined  b* 


(B-13) 


max 


( A ■ A ) 


'1 


and  .„(3)  denotes  the  maxi.mum  eigenvalue  of  B.  Applying 

Hi  3. 

these  two  properties  to  Eq . B-11  leads  to  the  result 


R' 


■{c 

I xz 


zz  xz ) 


(B-15) 


Combining 
''92..  ^ 


B-6,  B-10 , and  B-15  with  the  equality 
the  final  bound  becomes  simply 


(B-16) 


*For  a general  discussion  of  matrix  norms  and  their  properties 
see  Ref.  20,  Chapter  1 and  Ref.  12,  Section  2. 
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The  error  bound  in  Eq . B-16  is  awkward  to  compute 

since  it  involves  the  square  root  and  the  eigenvalues  of 

^ zz 

R through  Eq . B-14 . For  a symmetric  matrix  A,  a simpler  relation 
can  be  written 

llAlU  = 

where  the  matrix  norm  is 

riAii,_  = 

and  a. . are  the  elements  of 
ij 

simplex  error  bound  is 

£ < l|PRp-^!i2  < i!PRP~‘il^  (3-19) 

where  P is  an  arbitrary  non-singular  matrix  which  may  be 
chosen  to  minimize  the  bound. 

A useful  choice  of  P in  Ec . 3-19  is  P = , 

zz 

where  D is  a diagonal  matrix,  since  it  leads  to 

£ < IjD'  RD"'^i!^  (B-20) 

which  requires  only  the  com.putation  of  the  residual  matrix 

ik 

R,  plus  row  and  column  scaling.  For  the  norm  in  Eq . 3-18 
it  is  desirable  to  "balance"  the  matrix  so  that  all  (absolute) 
coltunn  sums  are  nearly  equal.  This  is  accomplished  approxi- 
mately by  defining  the  diagonal  matrix  D to  consist  of  the 
diagonal  elements  of  C.. 

y 


max 


(A)|  1 liAli^ 


(B-17) 


defined  by 
n 

max 

J i=i 


max  / la.. 


(B-18) 


A (Ref.  20).  It  follows  that  a 


‘Scaling  to  reduce  matri.x  norms  is  discussed  in  Ref.  12, 
Section  11. 
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For  illustration,  suppose  that  C is  obtained  from 
C by  striking  out  the  off-diagonal  elements.  (In  the  fast 
estimation  algorithm  C is  constructed  to  be  band  diagonal.) 
Then  since  C is  a covariance  matrix 

Z2 


D-  R D“^ 


= II  I - CT^C  IL 

1 " £ 2Z  0 " 1 


(ifj) 


(3-21) 


i=l 


where  the  r. . are  correlation  coefficients.  Thus  the  column 
sums  in  Eq.  B-21  are  all  of  the  same  order  of  magnitude  and 
bounded  by  (n-1). 


Interestingly,  the  error  bounds  in  Eqs . 3-16  and  5-20 
are  the  same  as  for  the  relative  error  in  a determ.inist ic 
linear  S3’stem.  Also  the  dependence  of  the  error  on  the  con- 
dition of  the  matrix  C.  is  easily  obtained.  For  example,  set 

^ 0 

P = in  Eq.  B-19  and  the  bound  can  be  written  (see  Eq . B-2) 


e 1 i!R|i2  = iic:^  oci'^  < <2(C,) 


!!5C|! 

ilCoil 


(3-22) 


where 


(C,)  = i!c 


-1 


'2^''c 


^ic  '! 

2 !'^o"2 


max 


(C.)/>  . (C.) 
' min  c 


is  the  condition  number  of  C-  in  the  euclidean  norm.. 
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